Compelling bounds on equilibration times—the issue with Fermi's golden rule

The last decades have seen a major progress in the field of equilibration in closed quantum systems [1]. Concepts like typicality [2–4] and the eigenstate thermalization hypothesis [5, 6] have been brought forth. Furthermore, it has been established that for an initial state ρ populating many energy levels, expectation values ⟨A(t)⟩ will generically be very close to their temporal averages for most times within the interval to which the average refers ('equilibration on average'). While the conditions for this statement to be true are rather mild concerning the observable A and the Hamiltonian H [7–9], the respective time interval may be very large. For specific observables it may, e.g. scale with the dimension of the relevant Hilbert space [10, 11]. Moreover, concrete examples are known in which the corresponding equilibration times for physically relevant observables scale as T_eq ∝ N^α, α ⩾ 1/2, where N is the size of the system. This result has been found for systems featuring long range [12] as well as short range interactions [13], albeit for a somewhat different definitions of equilibration times. In fact, already for mesoscopic many-body systems with standard interaction strengths, the required equilibration interval may be on the order of the age of the Universe [12]. Thus, although the above statements in some sense establish equilibration under moderate conditions in the very long run, it is unclear whether or not this equilibration will ever occur in a physically relevant period of time. Hence, the question of an upper bound on this relaxation timescale has recently been much discussed. Since it is always possible to find mathematically well defined, permissible initial states that fully exhaust the above, unsatisfactorily large time interval, most contributions focus on additional, physically plausible conditions. These conditions, which are intended to capture the actual, physical state of affairs, may be imposed on the initial state, the observable, the structure of the system, or combinations thereof [14–19]. In the present paper we primarily discuss results from [20]. The latter rests on assumptions on all of the above.

This paper is organized as follows: in section 2 we briefly present a main result from reference [20]. Furthermore, we elaborate on the lack of predictive power of this result in cases in which Fermi's golden rule applies. In section 3 we explain some models, each of which comprises a single spin in a magnetic field interacting with a (finite) bath, consisting of spins itself. We initialize the system in a standard system-bath product state and numerically solve the Schrödinger equation, monitoring the system's spin component parallel to the magnetic field. These data unveil the regime of validity of Fermi's golden rule with respect to the crucial system parameters. Section 4 discusses the scaling of the critical interaction strength at which open system predictions start to become unreliable. In section 5 the implications of the numerical findings from section 3 on the assumptions and statements from [20] are named and explained. Eventually, we sum up and conclude in section 6.

To begin with, we state the main result of [20] (hereafter called the García-Pintos bound (GPB)) in a comprehensive form. The GPB addresses an equilibration time T_eq. To further specify T_eq we introduce some notation. Let ρ be the initial state of the system. Furthermore, let A(t) denote an observable A in the Heisenberg picture and ⟨A(t)⟩ := Tr{A(t)ρ} its time dependent expectation value. Because the closed system dynamics is unitary (and the system is finite), ⟨A(t)⟩ has a well defined 'infinite time average' $\overline{\langle A\rangle }=:{\langle A\rangle }_{\text{eq}}$, which is routinely considered as the equilibrium value of A in case the observable A equilibrates at all [9]. Now, consider a deviation D(t) of the actual expectation value from its equilibrium, i.e. $D\left(t\right){:=}{\left(\langle A\left(t\right)\rangle -{\langle A\rangle }_{\text{eq}}\right)}^{2}/4{\Vert}A{{\Vert}}^{2}$, where ||A|| is the largest absolute eigenvalue of A. Moreover, consider an average of D(t) over the time interval [0, T] denoted by ${\overline{D}}_{T}$. The condition that defines T_eq is that ${\overline{D}}_{T}\ll 1$ must hold for T ≫ T_eq (for non-equilibrating systems such a T_eq may not exist [9]). The GPB is an explicit expression for such a T_eq (see equation (4)), based on ρ, A(0) and H, where H is the Hamiltonian of the system. As the GPB involves somewhat refined functions of the above three operators, we need to specify these before stating the GPB explicitly. A central role is taken by the probability distribution p_jk, which is defined as

$$\begin{equation}\begin{aligned}\hfill {p}_{jk}\propto \vert {\rho }_{jk}{A}_{kj}\vert \quad \text{for}\enspace {E}_{j}-{E}_{k}\ne 0,\\ \hfill {p}_{jk}=0\quad \text{for}\enspace {E}_{j}-{E}_{k}=0,\enspace \sum _{j,k}{p}_{jk}=1,\end{aligned}\end{equation} \tag{ 1 }$$

where E_j, E_k are energy eigenvalues corresponding to energy eigenstates |j⟩, |k⟩. Furthermore, matrix elements are abbreviated as ρ_jk := ⟨j|ρ|k⟩, A_jk := ⟨j|A(0)|k⟩. While the GPB is not limited to this case, here, we focus on the p_jk that may be described in terms of a probability density function w(G). All examples we present below conform with such a description and it is plausible that this applies to many generic many-body scenarios. Prior to defining w(G), we define w(G, ) as

$$\begin{equation}w\left(G,{\epsilon}\right){:=}\frac{1}{{\epsilon}}\sum _{j,k}{\Theta} \left(\frac{{\epsilon}}{2}-\vert {E}_{j}-{E}_{k}-G\vert \right){p}_{jk},\end{equation} \tag{ 2 }$$

where Θ is the Heaviside function. This definition is the standard construction of a histogram in which the p_jk are sorted according to their respective energy differences E_j − E_k. It is now assumed that there exists a range of (small but not too small) such that w(G, ) is essentially independent of variations of within this range. The w(G, ) from this 'independence regime' are simply abbreviated as w(G). Let the standard deviation of w(G) be denoted by σ_G. Additionally, let w_max denote the maximum of w(G). The quantities a and Q that eventually enter the GPB are now defined as

$$\begin{equation}a{:=}{w}_{\mathrm{max}}{\sigma }_{G},\quad Q{:=}\sum _{i,j:{E}_{i}\ne {E}_{j}}\frac{\vert {\rho }_{ij}{A}_{ji}\vert }{{\Vert}A{\Vert}}.\end{equation} \tag{ 3 }$$

We can now state the GPB:

$$\begin{equation}{T}_{\text{eq}}=\frac{\pi a{\Vert}A{{\Vert}}^{1/2}{Q}^{5/2}}{\sqrt{\vert \mathrm{Tr}\left(\left[\left[\rho ,H\right],H\right]A\right)\vert }}=\frac{\pi a{\Vert}A{{\Vert}}^{1/2}{Q}^{5/2}}{\sqrt{\vert \frac{{\mathrm{d}}^{2}}{\mathrm{d}{t}^{2}}\langle A\left(t\right)\rangle {\vert }_{t=0}\vert }}.\end{equation} \tag{ 4 }$$

Obviously, the GPB links T_eq to the initial 'curvature' of the observable dynamics ${\partial }_{t}^{2}\langle A\left(t\right)\rangle {\vert }_{t=0}$ (which is practically accessible, cf figure 7). An actual, concrete bound on the equilibration time by means of T_eq, however, only arises from equation (4) if the numerator can be shown to be in an adequate sense small or at least bounded. This is a pivotal feature on which the 'predictive power' of the GPB hinges. The crucial quantities in the numerator are a and Q. As it is practically impossible to calculate a from its definition for many-body quantum systems, García-Pintos et al instead offer an assumption.

They argue that a ∼ 1 may be expected for w(G) that are 'unimodal'. Unimodal means that w(G) essentially consists of one central elevation like a Gaussian or a box distribution, etc. Indeed, a is invariant with respect to a rescaling as w(G) → sw(sG), as it would result from rescaling the Hamiltonian as H → sH (here s is some real, positive number). García-Pintos et al also offer various upper bounds on Q for different situations.

In the remainder of this section, we explain in which sense the conclusiveness of the GPB is in conflict with Fermi's golden rule (FGR). Let us stress that this conflict does not concern the validity or correctness of equation (4) as such, the latter is undisputed. It only concerns the assumptions on a and Q, which are required to find an actual value or estimate for T_eq. (Note that there is some evidence (cf section 5) that specifically the assumption on a is violated, rather than the assumption on Q.) Consider a Hamiltonian consisting of an unperturbed part H₀ and a perturbation H_int.

$$\begin{equation}H={H}_{0}+\lambda {H}_{\text{int}}\end{equation} \tag{ 5 }$$

Furthermore, consider an observable A, which is conserved under H₀, i.e. [A, H₀] = 0. If H₀ has a sufficiently wide and dense spectrum and λ is small, FGR may apply under well investigated conditions [21–23]. The applicability of the FGR approach yields, in the simplest case, a monoexponential decay, i.e.

$$\begin{equation}\langle A\left(t\right)\rangle =\left(\langle A\left(0\right)\rangle -{\langle A\rangle }_{\text{eq}}\right){\text{e}}^{-t/{\tau }_{\text{rel}}}+{\langle A\rangle }_{\text{eq}},\end{equation} \tag{ 6 }$$

where τ_rel := rλ⁻² and r is a real, positive number that depends on H₀ and H_int. More refined approaches, such as the Weisskopf–Wigner theory or open quantum system approaches, also arrive at such exponential decay dynamics [24, 25]. In the relevant case ∂_t⟨A(t)⟩|_t=0 = 0, obviously equation (6) cannot apply at t = 0. In this case equation (6) is meant to apply after a short 'Zeno time' τ_zeno that is often very short compared to the relaxation time τ_rel [23]. (Note, however, that the denominator of equation (4) addresses a time below the Zeno time, if the latter is nonzero.) We now aim to find the principal dependence of quantities in equation (4) on the interaction strength λ. While the definition of T_eq as given at the beginning of the present section does not fix the relation of τ_rel and T_eq rigorously, for exponential decays it appears plausible to require at least

$$\begin{equation}{T}_{\text{eq}}{\geqslant}{\tau }_{\text{rel}}.\end{equation} \tag{ 7 }$$

For the denominator of equation (4) we find with equation (5)

$$\begin{equation}\sqrt{\left\vert \frac{{\mathrm{d}}^{2}}{\mathrm{d}{t}^{2}},\langle A\left(t\right)\rangle ,{\vert }_{t=0}\right\vert }=\sqrt{\vert {c}_{1}\lambda +{c}_{2}{\lambda }^{2}\vert },\end{equation} \tag{ 8 }$$

where c₁ = Tr([H_int, A][ρ, H₀]), c₂ = Tr([H_int, A][ρ, H_int]). Plugging equations (6)–(8) into equation (4) yields

$$\begin{equation}\pi a{\Vert}A{{\Vert}}^{1/2}{Q}^{5/2}{\geqslant}r\frac{\sqrt{\vert {c}_{1}\lambda +{c}_{2}{\lambda }^{2}\vert }}{{\lambda }^{2}}\end{equation} \tag{ 9 }$$

for the numerator of equation (4). Obviously, the numerator of equation (4) diverges in the limit of weak interactions, i.e. λ → 0. The latter holds even if c₁ = 0. This contradicts the central assumption behind the GPB as outlined below equation (4). Hence, the validity of FGR in the weak coupling limit and a conclusive applicability of the GPB are mutually exclusive. This is the first main result of the present paper. Note that, in the above scenario, the exponential decay is due to some small 'conservation breaking' part of the Hamiltonian. However, exponential decay also often occurs without some part of the Hamiltonian being particularly small, e.g. slow Fourier modes in sufficiently large systems. Although the practical success of FGR is beyond any doubt, the theoretical applicability of FGR rests on various assumptions on the system in question, so does the applicability of standard open system methods. In order to learn about the applicability of either the GPB or FGR from considering examples, we analyze some spin systems in the following section 3 by numerically solving the respective Schrödinger equations. This analysis is comparable to numerical investigations performed in [20]. However, other than García-Pintos et al we analyze the weak coupling limit and consider system sizes that are too large to allow for numerically exact diagonalization of the respective Hamiltonians.

While we analyze a number of concretely specified models below, it is important to note that these models just represent some generic instances of the system-bath scenarios which are routinely considered in open quantum system theory. The (non-integrable) baths share some properties with standard solid state systems, like periodicity and locality (in this respect they differ from the otherwise comparable models addressed in [26–28]). Other than that, the details of our modeling are not peculiar at all. We varied details of the bath Hamiltonians in piecemeal fashion and found all below results unaltered (cf appendix A.3). Our archetypal model is an isotropic spin-1/2 Heisenberg system, consisting of a single system spin coupled to a bath. The bath is rectangularly shaped with 3 × L spins and features periodic boundary conditions in the longitudinal direction resulting in a wheel-like structure (cf figure 1). Thus, the total number of spins is given by N = 3L + 1. The single system spin is subject to an external magnetic field in the z-direction and interacts with three neighboring bath spins in the transverse direction. This model is non-integrable in the sense of the Bethe-ansatz. The bath Hamiltonian reads

$$\begin{align}\hfill {H}_{\text{bath}}=& J\sum _{r=1}^{3}\sum _{i=1}^{L}\left({S}_{i,r}^{x}{S}_{i+1,r}^{x}+{S}_{i,r}^{y}{S}_{i+1,r}^{y}+{S}_{i,r}^{z}{S}_{i+1,r}^{z}\right)\hfill \\ \hfill & +J\sum _{i=1}^{L}\left({S}_{i,1}^{x}{S}_{i,2}^{x}+{S}_{i,1}^{y}{S}_{i,2}^{y}+{S}_{i,1}^{z}{S}_{i,2}^{z}\right)\hfill \\ \hfill & +J\sum _{i=1}^{L}\left({S}_{i,2}^{x}{S}_{i,3}^{x}+{S}_{i,2}^{y}{S}_{i,3}^{y}+{S}_{i,2}^{z}{S}_{i,3}^{z}\right),\hfill \end{align} \tag{ 10 }$$

where ${S}_{i,r}^{x,y,z}$ are spin-1/2 operators at site (i, r) and L + 1 ≡ 1. The exchange coupling constant J as well as ℏ are set to unity.

Zoom In Zoom Out Reset image size

The Hamiltonian of the system is given by

$$\begin{equation}{H}_{\text{sys}}=B{S}_{\text{sys}}^{z},\end{equation} \tag{ 11 }$$

where ${S}_{\text{sys}}^{x,y,z}$ denote the spin-1/2 operators of the additional system spin and B = 0.5. The interaction between bath and system is described by the Hamiltonian

$$\begin{equation}{H}_{\text{int}}=\left[\left({S}_{1,1}^{x}+{S}_{1,2}^{x}+{S}_{1,3}^{x}\right){S}_{\text{sys}}^{x}\right.\left.+\left({S}_{1,1}^{y}+{S}_{1,2}^{y}+{S}_{1,3}^{y}\right){S}_{\text{sys}}^{y}\right.\left.+\left({S}_{1,1}^{z}+{S}_{1,2}^{z}+{S}_{1,3}^{z}\right){S}_{\text{sys}}^{z}\right]\end{equation} \tag{ 12 }$$

and contributes with a factor λ to the total Hamiltonian

$$\begin{equation}H={H}_{\text{sys}}+{H}_{\text{bath}}+\lambda {H}_{\text{int}}.\end{equation} \tag{ 13 }$$

The considered initial states are product states of a system state π_↑ and a bath state π_E,δ. This corresponds to a situation where system and bath are initially uncorrelated and then brought into contact via H_int at t = 0. The system state is a projector onto the ${S}_{\text{sys}}^{z}$-eigenstate corresponding to spin-up. The bath state π_E,δ is a projector onto a (small) energy window of width δ centered around a mean energy E.

$$\begin{equation}\rho =\frac{{\pi }_{{\uparrow}}\otimes {\pi }_{E,\delta }}{\mathrm{Tr}\left\{{\pi }_{{\uparrow}}\otimes {\pi }_{E,\delta }\right\}}\end{equation} \tag{ 14 }$$

Concretely, we fix the width of the energy window δ = 0.1, which is very small compared to the scale of the full energy spectrum of the bath. Given the size of the systems, it comprises nevertheless a very large number of energy eigenstates. To keep track of finite size effects, we increment the baths circumference L in steps of size one, thus adding three spins to the bath in each step. An inverse temperature β is defined as ∂_E  log  Ω(E), where Ω(E) is the density of states of the bath at energy E. This 'microcanonical' definition of temperature is also employed in [20]. For comparability of different bath sizes we aim at keeping β fixed while incrementing the bath size. As H_bath is local, the bath energy is expected to scale linearly with the bath size. Hence, we choose a scaling of the initial bath energy as E ≈ −0.15(N − 1), which corresponds to choosing β ≈ 0.4. Given these specifications of the Hamiltonian and the initial state, we numerically solve the corresponding Schrödinger equation and monitor the expectation value of the z-component of the magnetization of the system-spin, i.e. $\langle {S}_{\text{sys}}^{z}\left(t\right)\rangle$. Some results are displayed in figure 2 for a schematic overview. In accord with open quantum system theory, these results suggest to distinguish three cases.

• (a)  
  Non-Markovian regime: for strong coupling (e.g. λ = 1.0) the magnetization quickly decays to the equilibrium value, i.e. ${\langle {S}_{\text{sys}}^{z}\rangle }_{\text{mc}}$, in a nonexponential way. The description of these dynamics requires the incorporation of memory effects in some way. While this is a very active field in open quantum theory, we do not investigate this regime any further in the present paper.
• (b)  
  Markovian regime: for weak coupling (e.g. λ = 0.1) there is a monoexponential decay to the thermal equilibrium. This exponential decay is in full accord with FGR. A large number of systems ranging from quantum optics to condensed matter fall into this regime [25, 29]. It also largely coincides with the field of quantum semi-groups and the Lindblad approach.
• (c)  
  Superweak coupling regime: for very weak coupling (e.g. λ = 0.01) the magnetization does not decay to the thermal equilibrium value at all, it rather gets stuck at a value closer to the initial value, which indicates the breakdown of FGR. This value depends on the interaction strength and on the bath size. In accord with standard open quantum system theory, our below results indicate that this regime only exists for finite baths. We are not aware of any systematic approach to this regime in the literature to date. As the conflict between the GPB and FGR arises in the limit of weak interactions, cf equation (9), we are primarily interested in the transition from the Markovian to the superweak regime. A prime indicator of superweak dynamics is, as mentioned, that $\langle {S}_{\text{sys}}^{z}\left(t\right)\rangle$ no longer decays down to the microcanonical expectation value ${\langle {S}_{\text{sys}}^{z}\rangle }_{\text{mc}}=-0.05$ as it does in the non-Markovian and the Markovian regime.

Zoom In Zoom Out Reset image size

Figure 3 shows the longtime average value of the magnetization plotted over the interaction strength λ for various bath sizes. For sufficiently strong coupling the magnetization decays to the thermal equilibrium value for all bath sizes. For each bath size there exists a critical interaction strength λ_crit, below which the magnetization gets stuck at a non-thermal longtime average value, thus signaling the transition from the Markovian to the superweak regime. This critical interaction strength λ_crit decreases with bath size. We chose $\overline{\langle {S}_{\text{sys}}^{z}\rangle }=-0.04$ (horizontal grey line) to define λ_crit. This choice is quite arbitrary, but it turns out that the below scaling of λ_crit is rather insensitive to the exact positioning of this threshold, as long as it is sufficiently close to the thermal equilibrium value. This finding is corroborated in appendix A.4. Obviously, one expects $\overline{\langle {S}_{\text{sys}}^{z}\rangle }\to 0.5$ as λ → 0. Figure 4 displays the critical interaction strength plotted over the inverse system size. It strongly suggests that λ_crit → 0 very quickly with increasing bath size N. Hence, for all mesoscopic to macroscopic systems, and even more so in the thermodynamic limit, the transition to the superweak regime practically never occurs, such that behavior other than Markovian can hardly be expected even for physically very weak interactions. While this finding is another main result of the quantitative analysis at hand, it qualitatively hardly comes as a surprise in a larger context, given the practical success of Markovian quantum master equations. However, to elaborate on this result somewhat further, we present a theory that captures the data in figure 4 rather accurately in section 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next we confirm the validity of FGR in the Markovian regime and discuss relaxation/equilibration times in all regimes. The motivation for the latter is twofold: on the one hand equilibration times enter the GPB (cf equation (4)), on the other hand the scaling of equilibration times with the interaction strength may serve as an additional, quantitative indicator for the validity of FGR. Figure 5 displays the observable dynamics $\langle {S}_{\text{sys}}^{z}\left(t\right)\rangle$ for all 11 interaction strengths and bath sizes from the Markovian regime, which is lower bounded by λ_crit as determined from figure 4 and upper bounded by λ_non-Mark ≈ 0.3 for all N. Note the time axis is scaled with the squared interaction strength such that a collapse of the data onto one decaying exponential indicates the accordance with equation (6) and hence FGR. This collapse is evident.

Zoom In Zoom Out Reset image size

In figure 6 the relaxation time τ_rel is plotted over λ⁻². Here τ_rel is the time at which the magnetization has decayed to 1/e of its original value relative to the equilibrium value (cf equation (6)). In the Markovian regime, i.e. for ${\lambda }_{\text{non}-\text{Mark}}^{-2}{\leqslant}{\lambda }^{-2}{\leqslant}{\lambda }_{\text{crit}}^{-2}$, the relaxation time scales as τ_rel ∼ λ⁻², as predicted by FGR, which also confirms the applicability of FGR in the Markovian regime. At very small λ, i.e. in the superweak coupling regime, τ_rel first increases more slowly with increasing λ⁻² and likely eventually even decreases to zero. From the data displayed in figure 6 this behavior is, however, qualitatively only visible for N = 16 due to numerical limitations at extremely small λ.

Zoom In Zoom Out Reset image size

Eventually, we directly numerically probe the connection between short time and long time dynamics suggested in equation (4). To this end we compute the 'initial curvatures' $\sqrt{\vert {\partial }_{t}^{2}\langle A\left(t\right)\rangle {\vert }_{t=0}\vert }$ for various interaction strengths and systems sizes. The result is displayed in figure 7. As expected from equation (8), and in full accord with a corresponding statement in [20], the square root of the initial curvature scales linearly with λ and is practically independent of the system size. We are now set to assess the crucial numerator from equation (4) numerically. From equations (4) and (7) follows

$$\begin{equation}\begin{aligned}\hfill \pi a{\Vert}A{{\Vert}}^{1/2}{Q}^{5/2}& ={T}_{\text{eq}}\sqrt{\left\vert \frac{{\mathrm{d}}^{2}}{\mathrm{d}{t}^{2}},\langle A\left(t\right)\rangle ,{\vert }_{t=0}\right\vert }\hfill \\ \hfill & {\geqslant}{\tau }_{\text{rel}}\sqrt{\left\vert \frac{{\mathrm{d}}^{2}}{\mathrm{d}{t}^{2}},\langle A\left(t\right)\rangle ,{\vert }_{t=0}\right\vert }.\hfill \end{aligned}\end{equation} \tag{ 15 }$$

The lower bound to the numerator is displayed in figure 8. Recall that for a conclusive application of the GPB this numerator must be appropriately upper bounded. Correspondingly, García-Pintos et al [20] offer estimates for both a and Q. While a ∼ 1 is simply traced back to the unimodality of w, the discussion on the order of magnitude of Q is quite involved. However, in the case of weak interactions, a microcanonical initial bath state comprising a large number of energy eigenstates, and an exponentially growing density of states in the bath (with an exponent β which is not too large), Q may also be expected to be of order unity, according to [20]. All these conditions apply to the models at hand. However, quite in contrast we find that the numerator grows at least up to values of ca. 55 already for N = 25 and interactions within the range of our numerical accessibility. Moreover, the data do not indicate any 'nearby' upper bound of the numerator at 55. This is at odds with a conclusive application of the GPB and another main result of the present paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Large numerators must occur, as outlined in section 2, for large systems in the Markovian regime on the verge to the superweak regime. Although figure 8 indicates that in the outer mathematical limit (which may be considered to be physically less relevant) λ → 0, the numerator may eventually be of order unity, its growth appears to continue substantially into the superweak regime. While we are unable to verify this directly, it appears plausible that the unbound growth of the numerator is due to an unbound growth of a, while Q ∼ 1 may very well hold. We argue for this plausibility in section 5.

We sum up this section as follows: in a potentially wide regime of interaction strengths, which is lower bounded by λ_crit, FGR is found to apply as suggested by standard open quantum system theory in the Markovian regime. The lower bound appears to decrease rapidly with system size N (cf also section 4), making it practically irrelevant for mesoscopic and macroscopic systems. This relates to the GPR in as much as the validity of FGR implies the breakdown of the practical applicability of the GPB at sufficiently weak interactions. We numerically confirmed the occurrence of this breakdown directly for a system comprising N = 25 spins.

While the results on λ_crit in figure 4 are model dependent, a similar scaling may be expected whenever H_int complies with the eigenstate thermalization hypothesis. This claim is substantiated in the following. The starting point is the assumption that within the superweak regime, the overlap of the eigenstates of the uncoupled system |n, sys + bath⟩ and those of the full system |n, sys + bath + int⟩ is relatively large, i.e.

$$\begin{equation}\vert \langle n,\text{sys}+\text{bath}\vert n,\text{sys}+\text{bath}+\text{int}\rangle \vert \approx 1.\end{equation} \tag{ 16 }$$

A strong indication for this to occur, arises from the leading order contributions to a perturbative correction to the eigenstates being small. This condition may, according to textbook level perturbation theory, be approximated as

$$\begin{equation}{\lambda }^{2}\sum _{m\ne n}\frac{\vert \langle m,\text{sys}+\text{bath}\vert {H}_{\text{int}}\vert n,\text{sys}+\text{bath}\rangle {\vert }^{2}}{{\left(\left(m-n\right)\frac{1}{{\Omega}\left(N,\beta \right)}\right)}^{2}}\ll 1,\end{equation} \tag{ 17 }$$

where Ω(N, β) is the density of states of a system comprising N spins (or other similar subsystems) at the energy that corresponds to the inverse temperature β. In equation (17) it is assumed that the level spacing within the relevant energy regime may be approximated as being constant. In this case the 'mean' level spacing is given by 1/Ω(N, β). Following the eigenstate thermalization hypothesis ansatz [30] we furthermore assume that there exists a typical value for the absolute squares of the matrix elements of the coupling operator, which varies with the energies E_n, E_m only on energy scales much larger than the one relevant here. Furthermore, the eigenstate thermalization hypothesis suggests a specific scaling of these matrix elements with the density of states. Following the eigenstate thermalization hypothesis we thus assume

$$\begin{equation}\vert \langle m,\text{sys}+\text{bath}\vert {H}_{\text{int}}\vert n,\text{sys}+\text{bath}\rangle {\vert }^{2}\approx \frac{{C}_{1}}{{\Omega}\left(N,\beta \right)},\end{equation} \tag{ 18 }$$

where C₁ is some real constant. Exploiting this, equation (17) turns into

$$\begin{equation}{\lambda }^{2}{C}_{1}{\Omega}\left(N,\beta \right)\sum _{k\ne 0}\frac{1}{{k}^{2}}\ll 1.\end{equation} \tag{ 19 }$$

As the sum assumes the finite value π²/3, we conclude

$$\begin{equation}{\lambda }_{\text{crit}}=\frac{{C}_{2}}{\sqrt{{\Omega}\left(N,\beta \right)}}\end{equation} \tag{ 20 }$$

for the scaling of λ_crit, where C₂ is a constant whose concrete value depends on H_sys, H_bath and H_int. Now we turn to an estimate for Ω(N, β). For a sufficiently large Heisenberg spin system (or any other system consisting of N similar, similarly and locally interacting subsystems) it is reasonable to assume that the density of states is Gaussian with mean zero and variance proportional to the particle number N.

$$\begin{equation}{\Omega}\left(N,E\right)\sim \frac{{2}^{N}}{\sqrt{N}}{\mathrm{e}}^{-{E}^{2}/\alpha N}\end{equation} \tag{ 21 }$$

Using β = ∂_E  log  Ω leads to

$$\begin{equation}{\Omega}\left(N,\beta \right)\sim \frac{1}{\sqrt{N}}{\mathrm{e}}^{\left(\mathrm{log}2-0.25\alpha {\beta }^{2}\right)N}.\end{equation} \tag{ 22 }$$

Inserting this into equation (20) and fitting for C₂ and α yields the dashed line in figure 4, which matches the data quite well. This remarkable agreement in turn backs up the argumentation which lead to equation (20). As the density of states is routinely expected to scale exponentially with the size of the system, equation (20) indicates that λ_crit will generally be exponentially small in the system size and thus the superweak regime will practically never be observed.

As explained in section 2 the applicability of the GPB hinges on the two parameters a and Q, both of which should be of order unity to establish a meaningful relation between the short time dynamics and the equilibration time in the sense of the GPB. However, for some of the numerical examples considered in section 3, at least one of the parameters must be substantially larger than unity. While we are unable to perform a direct numerical check for large system sizes, we strongly suspect that a ∼ 1 is violated at weak interactions even though the corresponding w(G) (cf section 2) is strictly unimodal. In the remainder of the present section we explain and back up this claim.

Consider the mathematically simple case of an infinite temperature environment in the initial state $\rho =\left({S}_{\text{sys}}^{z}+{1}_{\text{sys}}/2\right)\otimes {1}_{\text{bath}}/{d}_{\text{bath}}$. This initial state yields

$$\begin{equation}\langle {S}_{\text{sys}}^{z}\left(t\right)\rangle =\mathrm{Tr}\left\{{S}_{\text{sys}}^{z}\left(t\right)\rho \right\}=\mathrm{Tr}\left\{{S}_{\text{sys}}^{z}\left(t\right){S}_{\text{sys}}^{z}\right\}\end{equation} \tag{ 23 }$$

for the dynamics of the observable $\langle {S}_{\text{sys}}^{z}\left(t\right)\rangle$, which may be rewritten as

$$\begin{equation}\langle {S}_{\text{sys}}^{z}\left(t\right)\rangle =\sum _{j,k}\vert \langle j\vert {S}_{\text{sys}}^{z}\vert k\rangle {\vert }^{2}{\mathrm{e}}^{\text{i}\left({E}_{j}-{E}_{k}\right)t}.\end{equation} \tag{ 24 }$$

Now, consider the distribution p_jk as defined in equation (1) for this initial state and choice of observable, i.e. $A={S}_{\text{sys}}^{z}$.

$$\begin{equation}{p}_{jk}\propto \vert {\langle j\vert {S}_{\text{sys}}^{z}\vert k\rangle }^{2}\vert =\vert \langle j\vert {S}_{\text{sys}}^{z}\vert k\rangle {\vert }^{2}\quad \text{for}\enspace j\ne k.\end{equation} \tag{ 25 }$$

For non-integrable systems in the sense of a Bethe ansatz, the eigenstate thermalization hypothesis may be expected to hold, yielding $\vert \langle j\vert {S}_{\text{sys}}^{z}\vert j\rangle {\vert }^{2}\approx 0$. Exploiting this case, the insertion of equation (25) into equation (24) yields

$$\begin{equation}\langle {S}_{\text{sys}}^{z}\left(t\right)\rangle \propto \sum _{j,k}{p}_{jk}{\mathrm{e}}^{\text{i}\left({E}_{j}-{E}_{k}\right)t}.\end{equation} \tag{ 26 }$$

To the extend that p_jk may indeed be replaced by a smooth probability density as discussed around equation (2), equation (26) may be rewritten as

$$\begin{equation}\langle {S}_{\text{sys}}^{z}\left(t\right)\rangle \propto \int w\left(G\right){\mathrm{e}}^{\text{i}Gt}\mathrm{d}G.\end{equation} \tag{ 27 }$$

Thus, for the present scenario, w(G) is essentially the Fourier transform of the observable dynamics $\langle {S}_{\text{sys}}^{z}\left(t\right)\rangle$. Based on the numerical findings displayed, e.g. in figure 2 it appears plausible that $\langle {S}_{\text{sys}}^{z}\left(t\right)\rangle$ will be an exponential decay for infinite temperature initial states as well. Therefore, w(G) will be Lorentzian. While a Lorentzian distribution is clearly unimodal with one well-behaved maximum, its variance diverges. Consequently a, as defined in equation (3), diverges as well. Thus, in contrast to the assumptions in [20], a ∼ 1 does not hold. This is the last main result of the present paper. Of course σ_G cannot really diverge in any system featuring a finite energy spectrum. However, the finiteness of the spectrum essentially causes a cut-off of the tails of the Lorentzian at some frequency. This cut-off actually renders the standard deviation σ_G finite. Nonetheless, this standard deviation does not reasonably reflect the width of w(G). It will be much larger than other measures of the width such as the full-width-at-half-maximum, etc.

Some attention should also be paid to the question whether or not the GPB scales with the size of the environment. (Earlier works presented upper bounds that explicitly depend on the size of the environment, which is often seen as a drawback [10, 11].) While the GPB does not explicitly dependent on the size of the environment, the latter may enter via the parameter a. For any (weak) interaction strength λ there exists a system size N(λ_crit) above which FGR applies. Above that size the GPB is independent of N. Below or at N(λ_crit), however, the numerator in equation (4) may depend on N rather strongly. For arbitrarily small λ this N(λ_crit) may become arbitrarily large. Thus, in the class of models discussed in the paper at hand, one can always find instances for which the GPB depends on system size even for very large systems, i.e. N ≫ 1.

In the paper at hand we conceptually and numerically analyzed an upper bound on equilibration times presented in a recent paper by García-Pintos et al To this end, we investigated a standard system-bath setup by monitoring the system's magnetization for various bath sizes and interaction strengths. This numerical investigation is based on the solution of the time dependent Schrödinger equation for the full system, including the bath. We identified a Markovian regime of interaction strengths λ in which Fermi's golden rule holds, i.e. the system thermalizes in an exponential way and the equilibration time scales as λ⁻². This relates to the García-Pintos bound inasmuch as the validity of Fermi's golden rule and the usefulness of the García-Pintos bound are analytically shown to be mutually exclusive at sufficiently small λ. At extremely small λ, we indeed find a 'superweak' regime in which Fermi's golden rule does not apply. This regime (in principle) exists for finite baths and is reached below some λ_crit which is shown to scale inversely exponentially in the bath size, suggesting that the superweak regime practically ceases to exist when considering moderately large systems. However, in the superweak regime the García-Pintos bound may eventually regain applicability, although its non-applicability is found to extend also into the superweak regime.

Stimulating discussions with R Steinigeweg and J Richter are gratefully acknowledged. We also thank L P García-Pintos and A M Alhambra for interesting discussions, as well as A J Short for a comment on an early version of this paper. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) within the Research Unit FOR 2692 under Grant No. 397107022.

Details of our numerical implementation will be discussed in this appendix. We make use of the concept of typicality (cf appendix A.1) and use a time evolution algorithm (real and imaginary) based on Chebyshev polynomials (cf appendix A.2). Since the full Hamiltonian conserves magnetization, we perform the procedure outlined below in each magnetization subspace. The dynamics in the full Hilbert space are obtained by piecing together the contributions of each magnetization subspace weighted with their respective binomial weight. Note that the main hindrance to our calculations is not the exponentially large Hilbert space dimension, but rather the extremely long times that have to be reached in real time for small interaction strengths. In appendix A.3 a result for randomized bath couplings is shown.

A.1. Typicality

As mentioned in section 3, the initial state is a product state of a microcanonical bath state and a projected spin-up system state. Since the numerical integration of the von-Neumann equation can be cumbersome, we make use of the concept of typicality, which states that a single 'typical' pure state can have the same thermodynamic properties as the full statistical ensemble [31]. Not only is it more memory efficient to work with pure states, the availability of efficient time evolution algorithms for pure states, e.g. Runge–Kutta or Chebyshev polynomials, constitutes a major advantage. To find such a typical state, a pure state |ϕ⟩ is drawn at random from the Hilbert space according to the unitary invariant Haar measure.

$$\begin{equation}\vert \phi \rangle =\sum _{i}{c}_{i}\vert i\rangle \end{equation} \tag{ A.1 }$$

Real and imaginary part of the complex coefficients c_i are drawn from a Gaussian distribution with mean zero and unit variance and the set {|i⟩} is an arbitrary basis of the Hilbert space, e.g. the Ising basis. Consider the new normalized state

$$\begin{equation}\vert \psi \rangle =\frac{\sqrt{\rho }\enspace \vert \phi \rangle }{\sqrt{\langle \phi \vert \rho \vert \phi \rangle }}.\end{equation} \tag{ A.2 }$$

It can be shown [32] that for the overwhelmingly majority of random states |ϕ⟩, the pure state |ψ⟩ exhibits effectively the same thermodynamic behavior as the mixed state ρ, i.e.

$$\begin{equation}\langle A\left(t\right)\rangle =\mathrm{Tr}\left\{\rho A\left(t\right)\right\}=\langle \psi \vert A\left(t\right)\vert \psi \rangle +\varepsilon .\end{equation} \tag{ A.3 }$$

Importantly, the induced error ɛ = ɛ(|ψ⟩) has mean zero, i.e. $\overline{\varepsilon }=0$, and a standard deviation that scales inversely proportional to the square root of the effective Hilbert space dimension, i.e. $\sigma \left(\varepsilon \right)\propto 1/\sqrt{{d}_{\text{eff}}}$ [33]. The effective dimension d_eff = 1/Tr{ρ²} is a measure of how many pure states contribute to the mixture ρ.

Zoom In Zoom Out Reset image size

In the paper at hand the initial state ρ (cf equation (14)) is a projection operator. Therefore, it is permissible to drop the square root in the numerator in equation (A.2). Now the projectors π_↑ ⊗ 1 and 1 ⊗ π_E,δ need to be applied to the state |ϕ⟩, which is an element of the product Hilbert space. As we are working in the Ising basis, the action of π_↑ ⊗ 1 on |ϕ⟩ is easily implemented by setting corresponding components of the state vector to zero. Since it is unfeasible to diagonalize the full many-body Hamiltonian, we replace the bath projector by a Gaussian filter which suppresses contributions of energy eigenstates far away from the desired energy E, resulting in a narrowly populated energy window of width (variance) δ.

$$\begin{equation}{\pi }_{E,\delta }\approx \mathrm{exp}\left(\frac{-{\left({H}_{\text{bath}}-E\right)}^{2}}{2\delta }\right)\end{equation} \tag{ A.4 }$$

The Gaussian filter is applied with a Chebyshev algorithm (cf appendix A.2). Lastly, the resulting wave function is normalized to obtain the state |ψ⟩ as in equation (A.2). In this scenario the effective dimension d_eff is essentially the number of states in the energy window. Increasing the size of the system, while keeping δ fixed, results in an exponentially growing effective dimension d_eff and therefore in a negligible typicality error for moderately sized systems.

A.2. Chebyshev polynomials

A Chebyshev type algorithm is employed in order to evolve a pure state in real and imaginary time [34, 35]. Say it is desirable to approximate a scalar function f(x) in the interval [−1, 1] by a polynomial expansion, i.e.

$$\begin{equation}f\left(x\right)\approx \sum _{n}{d}_{n}{P}_{n}\left(x\right),\end{equation} \tag{ A.5 }$$

with coefficients d_n and polynomials P_n of order n. The unique set of polynomials that minimizes the maximum error in this interval is called Chebyshev polynomials of the first kind. They are denoted by T_n(x) and can be written down recursively as

$$\begin{equation}{T}_{n}\left(x\right)=2x{T}_{n-1}\left(x\right)-{T}_{n-2}\left(x\right)\end{equation} \tag{ A.6 }$$

with T₀(x) = 1 and T₁(x) = x. They are orthogonal with respect to the weighted scalar product

$$\begin{equation}\langle {T}_{n}\vert {T}_{m}\rangle {=\int }_{-1}^{1}\frac{{T}_{n}\left(x\right){T}_{m}\left(x\right)}{\pi \sqrt{1-{x}^{2}}}\enspace \mathrm{d}x={\delta }_{nm}{D}_{n}\end{equation} \tag{ A.7 }$$

with D₀ = 1 and D_n>0 = 1/2. In order to approximate the time evolution operator, the bandwidth of the Hamiltonian has to be rescaled accordingly. Defining ω = (E_max − E_min)/2 and χ = (E_max + E_min)/2, where E_max (E_min) is the maximal (minimal) energy eigenvalue, the rescaled Hamiltonian is obtained by $\tilde {H}=\left(H-\chi \right)/\omega$. In practice a small safety parameter is chosen that ensures that the rescaled spectrum lies well within [−1, 1]. Now we get

$$\begin{equation}{\mathrm{e}}^{-\mathrm{i}H{\Delta}t}={\mathrm{e}}^{-\mathrm{i}\chi {\Delta}t}\left[\right.{d}_{0}\left(\omega {\Delta}t\right)+2\sum _{n{\geqslant}1}{d}_{n}\left(\omega {\Delta}t\right){T}_{n}\left(\tilde {H}\right)\left.\right]\end{equation} \tag{ A.8 }$$

with complex coefficients

$$\begin{equation}{d}_{n}\left(\omega {\Delta}t\right){=\int }_{-1}^{1}\frac{{T}_{n}\left(x\right){\mathrm{e}}^{-\mathrm{i}x\omega {\Delta}t}}{\pi \sqrt{1-{x}^{2}}}\mathrm{d}x={\left(-\mathrm{i}\right)}^{n}{J}_{n}\left(\omega {\Delta}t\right),\end{equation} \tag{ A.9 }$$

where J_n denotes the nth order Bessel function of the first kind. Since the coefficients only depend on the time step, but not on time itself, they only have to calculated once. Applying equation (A.8) to a state |ψ(t)⟩ boils down to calculating ${T}_{n}\left(\tilde {H}\right)\vert \psi \left(t\right)\rangle$ for various n, which can be done iteratively using equation (A.6). Terminating the sum in equation (A.8) at an upper bound M gives the Mth order Chebyshev approximation of the time evolution operator. The required order for convergence depends on the particular problem. For our biggest system with N = 25 spins we had to go up to order 40.

A.3. Validity of numerical results for a larger class of systems

A.4. Scaling of fit parameters with choice of critical threshold

The choice $\overline{\langle {S}_{\text{sys}}^{z}\rangle }=-0.04$ to define λ_crit is rather arbitrary. However, figures A.2 and A.3 show that the fit parameters C₂ and b are quite insensitive to the exact positioning of this threshold, as long as it is sufficiently close to the thermal equilibrium value. Furthermore, the principled form of the curve in figure 4 is quite robust to small variations of the fit parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size